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systems. This allows us to find the explicit formulas for all fitting parameters 
such as decay rates, coupling coefficients and characteristic intensities in terms 
of the mode profiles. The advantages of using the semi-analytical approach 
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1. Introduction 

The phenomenological Coupled Mode Theory (CMT) has been widely applied in optoelec- 
tronics, photonics and quantum optics to describe the linear and nonlinear properties of 
resonators [l}|3]. It played an importation role in the development of devices which can be 
embedded into photonic crystals: from waveguide splitters and add-drop filters [4|[5] to op- 
tical diodes and transistors (6|[7]. Moreover, it was used to achieve an efficient generation 
of harmonics and difference frequencies in microcavities [8l|9|, to describe bistable and mul- 



tistable switching 10,11 , and to explain self-pulsations and chaotic behavior in coupled 
microcavities 12 , 13 . The primary advantage of CMT comes from the fact that it allows one 
to understand properly the physical interactions between different modes of a system which 
are often hidden in the strictly numerical simulations. 
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The CMT equations can be obtained from general physical concepts like conservation 



of energy and time-reversal symmetry 14,15 . As a consequence, these equations contain 



several fitting parameters: decay rates of resonances, coupling coefficients to different scatte- 
ring channels and characteristic powers which describe the strength of the nonlinear effects. 
These parameters can be extracted from the experimental data, but most often they are de- 



termined after performing additional simulations in the time domain 16 17 . The main goal 
of this paper is to show that the CMT equations for the three-dimensional microcavities can 
be derived directly from the Maxwell equations without resorting to the phenomenological 
concepts. As a result, we were able to obtain the explicit formulas for all fitting parameters 
in terms of the mode profiles. We restrict our attention to the on-channel (or resonantly 
coupled) microcavities as opposed to the off-channel (side-coupled) microcavities 18 , since 
the former case can be readily applied to multilayered structures, and the results can be 
simplified considerably. 

The paper is organized as follows. In Section [2| we formulate the eigenvalue problem for 
a microcavity with two coupling ports and construct a complete set of orthogonal modes to 
expand an arbitrary field in it. To treat the electric and magnetic fields on equal footing, 
the Maxwell equations are written in a form which is similar to the Schrodinger equation. 
It is emphasized that due to the time-reversal symmetry the Maxwell equations are doubly 
degenerate, and two fundamental modes exist for any resonant frequency. In Section [3} we 
use these modes as a basis to describe the behavior of the microcavity in the vicinity of 
resonance. The CMT equations are derived for both cases when the fundamental modes 
are represented by traveling and standing waves. It is shown that the standing waves basis 
is particularly suitable for microcavities of high quality factors. Section |4] explains how to 
take into account perturbations caused by the Kerr nonlinearity and how to define the 
transfer matrix for the nonlinear microcavities. It provides the full set of equations for the 
time domain and frequency domain simulations including the explicit formulas for all fitting 
parameters. Section |5] presents several numerical examples and compares the accuracy of the 
CMT equations with other methods. It is also demonstrated how to apply the CMT equations 
to describe the nonlinear properties of microcavities with several localization centers. 

2. Generalized eigenvalue problem for Maxwell's equations 

2. A. Analogy with the Schrodinger equation 

To work with the electric and magnetic fields on equal footing, the Maxwell equations 

cV X E = -MH, (1) 
cV X H = edtE, (2) 
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where e is permittivity, /i is permeability, and c is the speed of hght in the vacuum, can be 
rewritten in a form which is similar to the Schrodinger equation 

ihpdt\^) = t\^) . (3) 

The wave function combines the electric and magnetic fields in a single vector 



I*) 



E 
H 



(4) 



(5) 



with the inner product defined as 

(vl/,|vl/,)= JJJ(E:.E, + H:-H,)dl-, 

where the integration is performed over the volume of the microcavity as shown in Fig. [l](a). 

The operator p is used as a weighting function, and it takes the following form for isotropic 
materials 

^ eI 


where I is the identity tensor. The generalization to the case of more complicated constitutive 



(6) 



relations such as those present in anisotropic or bi- anisotropic media is straightforward 19 
. It is important that the weighting operator is Hermitian (p^ = p) for lossless media. 
The operator L can be interpreted as a Hamiltonian and is defined as 
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ihc 



Vx 
-Vx 



(7) 



To prove that this operator is Hermitian, it is sufficient to show that (\l/a|L^6) = ijj'^ a\^b) for 
two arbitrary solutions with the same boundary conditions. This relation can be transformed 
to a surface integral around the microcavity [Fig. [l|^a)] by using the identity div(A x B) = 
B rot A - A rot B 

jJ(E:xH, + E, xH:)ds = 0. (8) 
The integral is nonzero only at the input and output ports so that it can be reduced to 



(e: X h, + e, X h:; 



\x=L 
\x=0 



0. 



(9) 



This equation is always valid for solutions that satisfy the periodic (Bloch) boundary condi- 
tions at the couphng ports 



\^!{x = L)) =e^'^|^(x = 0)) 



(10) 



where is a real constant. 
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2.B. Orthogonality relations between modes 

Assuming the time dependence of the form |\&(r, t)) = \ip{r)) e~"^*, the generahzed eigenvalue 
problem for Eq. (|3| can be formulated as 

L IV'm) = hUraPli^m) ■ (H) 

Since both operators L and p are Hermitian, the eigenvalues, or resonant frequencies Um, 
should be reaL Moreover, the eigenvectors, or modes corresponding to different frequencies 
Un 7^ oJmi should Satisfy the orthogonahty relation 

{ipnVP'ipm) = JJJ (^E: ■ E„ + /iH; ■ U^)dV = 0. (12) 

These modes can be normalized to a volume of the microcavity, but we prefer to preserve their 
meaning as the total energy stored in the mode and will normalize the boundary conditions 
instead. 

It is useful to introduce another Hermitian operator 



/ 
-/ 



(13) 

which satisfies the following commutation relations ct^L = — Lct^ and &zP = p&z- By using 
these properties, it can be shown that the solutions of Eq. (11) come in pairs, and for any 
solution with a positive frequency one can immediately construct another solution with a 
negative frequency 

L \a-z1pm) = -hUmP \o-ziJm) ■ (14) 

It is possible to return back to the positive frequencies by applying the operation of complex 
conjugation and to obtain two independent solutions for the same frequency Um 

m = ( hI ) = ( -nt ) ■ ^''^ 

In homogeneous media, these solutions can be considered as plane waves moving in the 
opposite directions {Ey, ±Hz)^ exp[i{±kxX — ut)]. 

Similar to the derivation of the orthogonality relations, it can be proved that for Un 7^ — 

M^.M = JJJ(^E; ■ - /iH: ■ UjdV = 0. (16) 

For modes of the same frequency cOn = cOm, it shows that the averaged electric energy stored 
in the mode is equal to the magnetic energy JJJ e|EmpdV^ = JJJ /i|Hmpdl^. If Eqs. (12) and 
(16) are combined, the orthogonality relations can be separately formulated for the electric 
and magnetic fields {un 7^ iwm) 

JJJ bE: . E^dV = JJJ /iH; ■ U^dV = 0. (17) 
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3. Slowly varying envelopes 

3. A. Traveling waves basis 

In the vicinity of the resonant frequency Um, the field can be searched in the following form 

\^) = a+(r,t) |^+>e-^--* + a_(r,t) \i^-)e-'--\ (18) 

where a±{r,t) are slowly varying envelopes of forward and backward moving waves. Substi- 
tuting the approximate solution (18) into the Schrodinger equation ^ gives 



ihp {dta+ |^+) + dta- |^^)) = L(a+) |^+) + L(a_) 



(19) 



where the notation L(a-i-) means that the derivatives in the operator L are applied only to 
the scalar envelopes a± 

(Va)x 
-(Va)x 

Projecting Eq. (19) on the vectors gives a set of two coupled equations which describe 
the propagation of the envelopes 



Ua) 



ihc 



(20) 



dta+ + (vg ■ V)a+ = -g*dta. 
dta^ — (vg ■ V)a- = —gdta^ 



(21) 
(22) 



The coupling terms appear due to the fact that the two fundamental modes 
orthogonal 

(V'-|«>=2jJJeE^d\/. (23) 
It is convenient to measure their overlap by using their common norm 



£|E„PdV 



and to introduce a special parameter 



(24) 



(25) 



Since the inequality J f{x)dx\ < J\f{x)\ dx holds for any complex function /(x), the 
parameter g is limited to the range \g\ < 1. It tends to zero in homogeneous media due 



to the rapid oscillation of phase in the numerator of Eq. (25). As a result, the coupling 



terms in Eqs. (21)-(22) disappear, and the two envelopes propagate independently. On the 



contrary, the field in nonhomogeneous media is localized around the defect regions, the 
amplitude varies very quickly in comparison to the phase, and this leads to large values of 
the parameter g. As a consequence, g can describe the localization strength of a resonance. 
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In the derivation of Eqs. (21)-(22), it was also used that 

= T2ihcVa ■ JJJ Re(E:;, x H„)dV^, 



0, 



(26) 
(27) 



and the parameter Vg was introduced 



JJJ Re(E:;,xHjd\/ 



(2^ 



The integrand in the numerator of Eq. (28) contains the averaged energy flow Siy 



(c/Svr) Re[E x H*], which is a solenoidal vector field div Sw = for any resonant frequency, 
because the averaged energy density W = (l/167r)(e|Ep + /i|Hp) does not change in a sta- 
tionary state. In multilayered structures, the energy flow depends only on the x-coordinate, 
which means that divSv^ = dx{Sw)x = 0, and as a result {Sw)x = const can be easily in- 
tegrated. Therefore, the parameter Vg can be interpreted as an effective group velocity. The 
same considerations should hold even for more complicated structures if the center of the 
microcavity is on the x-axis. Due to the symmetry, the dominant energy flow is also directed 
along the x-axis, and the directional derivative in Eqs. (23)-(24) can be approximated as 
(vg ■ V) = (vgd^). 

It is convenient to choose the position of reference planes in such a way that the boundary 
conditions (10) have = 0. Using these reference planes as channels for the in- and outgoing 



waves, it is possible to develop a scattering formalism [Fig. [T](b)] and to find the transmission 
spectrum of the microcavity in the vicinity of the resonance Um- Assuming that the reference 
planes are located at x = and x = L, the boundary conditions for the forward and backward 
moving envelopes can be written as a±(0) = u±, a±{L) = f=p. Approximating the spatial 
and temporal derivatives in Eqs. (21 )-([22)) with finite differences dxa± = (f=p — u±)/L and 
dta± = —iSu){u± + v^)/2, one can show that the scattering matrix Suv defined as 





ru t 
t r„ 




(29) 



up to the first order of the detuning from the resonance 6uj = u — Um is 

1 



where 7 = fg/L or 



1 — i{6uj/'-)) 
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ig{6uj/-f) 1 
1 i9*{Suj/-f) 



ca 



(30) 



(31) 
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which is obtained from Eq. (28) in the assumption that the mode profiles are considered as 
dimensionless quantities and normahzed in such a way that 

a = JJ Re[(E^:; x llj^]dydz = 1 cml (32) 

The transmission t{u) and reflection r^^^lu) spectra for waves incident on the ports U 
or V can be found by direct comparison of the matrix elements in Eqs. (29) and (30). For 
example, the transmission coefficient is given by 

t{co) = [l-ti5u/^)]-\ (33) 

and thus the resonance contour has the Lorentzian shape with the half-width at half- 
maximum equal to 7. 

3.B. Standing waves basis 

As was mentioned before, the basis formed by two modes {ipm) is not orthogonal. However, 
a linear combination of these modes can be used to construct a new basis which will have 
such a property. This is particularly easy to do for mirror symmetric structures 

(34) 
(35) 

The new modes {iPa) and \iPb) correspond to standing waves because the energy flow Sw for 
them is zero by definition. It can be shown that when one of them is exponentially growing 
around a defect region, the other mode is exponentially decaying [Fig. [2]. As a result, the 
norms of these modes can differ significantly for resonances with strong localization. By using 
the property lip^) = IV'a) ± ^ IV'b), the parameter g in Eq. (25) can be rewritten as 

_ JJJ.(Ei-E|)dl^ 
' JJj£(Ei + E|)dV'- 

Therefore, g is real for mirror symmetric structures and tends to ±1 depending on which 
mode dominates. 

The total field can be expanded in the new basis as 

1^) = A(r, t) ItfjA) e"^'^'"* + tB{r, t) {i/jb) e"*"™*, (37) 

where A{r,t) and B{r,t) are slowly varying amplitudes of the modes. The equations which 
describe their propagation can be found by substituting Eq. ( [37| ) into Eq. ^ 

ihp {dtA \^a) + tdtB I^b)) = HA) I^a) + tHB) I^-b) • (38) 
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Making projection of this equation on {iPaI and then on {ip^l leads to 



{l + g)dtA + {Y^-V)B = 0, 
(1 - g)dtB + (vg ■ V)A = 0, 

where the following properties were used 



A 



L(A)7Aa) = (V^B L(5)V^b) =0 



,(5)^b) = -hcVB . JJJ Re(E:; x UjdV. 



(39) 
(40) 



(41) 
(42) 
(43) 



The Eqs. (39)-(40) can be applied to develop the scattering formalism in the standing 
waves basis. The boundary conditions can be obtained by using the following relations with 
the traveling waves basis: A = a+ + a_ and i? = a+ — a_ . It is convenient to consider A and 
-B as a function of time only by taking their average value on the boundaries 



A 



B 



B{0) + B{L) 



w 



hich can be written in the matrix from as 



1 

1 



A 



The spatial derivatives can be approximated then by 

dA _ A{L)-A{0) _ ^ B- 
dx L 



B 



(44) 
(45) 

(46) 



L 

OB _ B{L) - 5(0) _^A- (n+ 
dx L L 



Therefore, Eqs. (39)-(40) take the following form 

l + gdA 



2 dt 
l~gdB 

2 dt 

where the decay rate 7 is given by 

7 



7A + 7(m+ + t;+), 
-fB + -f{u+ -v+). 



ca 



(47) 
(48) 

(49) 
(50) 

(51) 



JIf.(Ei + E|)dl-- 

A considerable simplification can be made for resonances with large quality factors. In this 
case, the parameter g tends to ±1 depending on which mode dominates. If the even mode 
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has a larger norm JJJ eE^dV^ ^ JJJ £:E|dV^, then according to Eq. (36) g ^ 1 and Eq. (50) 
reduces to i? = n+ — v^. The full set of the CMT equations is thus 



dA 
~dt 



U- 



1 



A 



1 



(52) 
(53) 



On the contrary, if the odd mode has a larger norm, then (7 — )■ — 1 and Eq. (49) reduces to 
A = + v+. The full set of the CMT equations is thus 



dB 
~dt 

V- 



-lB + -f{u+-v+), 



v+l \ 1 



(54) 
(55) 



4. Perturbations caused by the Kerr nonlinearity 

4- A. Time domain 

The nonlinear effects can be treated as perturbation terms in Eq. ([3]), and the standing 
waves basis is particularly suitable for this purpose because the dynamics of the system can 



be described only by one variable. The left hand side of Eq. (38) should be modified to 
include the perturbation term 



huj^Ap{A\i)A) +iB\i)B)) ■ 



(56) 



As a source of Ap, the Kerr nonlinearity will be considered, which corresponds to the fol- 
lowing constitutive relation between the electric field and the displacement vector 



D = eE + £kE^ 



(57) 



If the effect of the third harmonic generation, which is represented by the first term in the 
expansion 

(Re[E<,e-'"*])3 = (1/4) ReiE^e'^^'^*] + (3/4)|E,|2 RefE^e"*"*], (58) 
can be neglected, the constitutive relation is reduced to 



D = eE + (3/4)eK|E|2E, 



which means that 



£K|EpJ 




(59) 



(60) 
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Since the electric field can be written as E = AEa + ^-BEb, the projection of Eq. (56) on 
(V'aI and then on (?/'b| leads to the overlap integrals of the following form 



APB9 



3a;^ JJJeK(EAr(EB)W 
4 JJJ 5(Ei + E|)dr ' 



(61) 



where p and q are nonnegative integers with an additional restriction p + q = 4. There is 
only a small number of nonzero coefficients F which should be taken into account due to the 



fact that Ea and Eb are functions of the opposite parity. Eqs. (39)-(40) in presence of the 
Kerr nonlinearity are consequently 



:i + g)dtA + (vg ■ V)B = tA\A\^T AAAA + i{2A\B\^ - A*B')Taabb, 
;i - g)dtB + (vg ■ V)A = iB\B\''Tbbbb + ^(25|A|2 - B*A^)Taabb. 



(62) 
(63) 



The nonlinear effects described by these equations include the self-phase modulation (terms 
proportional to and -B|-Bp) as well as the cross-phase modulation (A|Sp and -B|y4p). 

If the coefficients F were equal, the relative strength of these effects would be given by the 



factor of 2, which coincides with the results obtained in the case of shallow gratings 21 . It is 
worth noting that there are a few additional terms {A*B^ and B*A'^), which are responsible 
for the phase conjugation and do not exist in the shallow gratings [22] . 

It turns out however that only the self-phase modulation is important for resonances with 
large quality factors. If the even mode \iPa) dominates, the overlap integrals differ significantly 
Taaaa ^ Taabb ^ Tbbbb- Therefore, to take into account the Kerr nonlinearity the CMT 



equations (52) and (54) should be modified in the following way 



dA 
~dt 
dB 
~dt 



7- 



7- 



Ia 
Bl_ 

Ib 



7 



7 



v4 + 7(m+ + 'y+), 
B + -i{u+-v+). 



(64) 
(65) 



The main infiuence of the Kerr nonlinearity results in the shift of the resonant frequency. To 



emphasize this, the time dependence exp(— iwrn^) used as a factor in Eq. (37) was explicitly 



included in the amplitudes A and B. Two new parameters were introduced I a and Jb which 
have the meaning of characteristic intensities. Before giving the explicit formulas for them, 
it is important to choose proper units for the electric field. 

The intensities will be measured in MW/cm^, and the nonlinear refractive index n2 caused 
by the Kerr nonlinearity will be specified in cm^/MW for consistency. The transition from the 
intensity dependent refractive index n{I) = n + to the permittivity can be performed 
as s{E) = e + 2en2\E\'^, where the units for the electromagnetic field were defined so as 
to produce a unit energy fiow in vacuum, namely /unit = (c/87r)|£'unitp = 1 MW/cm^. 
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Comparing with Eq. (59) gives £k = (8/3)en2 Therefore, the characteristic intensities are 

ccr 



'A,B 



(66) 



4-3. Frequency domain 



It is worth noting that the CMT equations (64) and (53) for resonances with even modes (or 



Eqs. (65) and (55) for odd modes) represent an extension of the scattering matrix method 



to the time domain. It is convenient to combine these equations 

Up- 



dA 



- M Wo - 7 





(67) 
(68) 



where A{t) denotes the amplitude of the dominating mode and the upper (lower) sign should 
be used if this mode is even (odd). The scattering matrix in the frequency domain is 



1 



i((5wcfr/7) 



±i{6uJcs/l) 1 

1 ±i{6ucs/l) 



(69) 



where the effective frequency detuning was introduced 6ues = uj — uq + 7|y4p//o which takes 
into account the shift of the resonant frequency due to the influence of the Kerr nonlinearity. 
The amplitude A depends on the amplitudes of the ingoing waves and can be found as a 
solution of the following equation 



1 



c^o ^ \A\ 



(70) 



It is a cubic equation which has three different roots in general case. This agrees with the 
fact that the system can show several stable states for the same input signals and explains its 
bistable behavior from the mathematical point of view. It is also possible to find A by using 



Eq. (68), which does not involve any nonlinear equations, however A becomes a function 



of both ingoing and outgoing signals. This can be particularly suitable when signals are 
incident only from one side. For example, if the incidence from the left is considered, f + = 



and A = — 1>_. Since the transmitted and input powers are given by Po 



a\V- 



and 



Pin = cr\u^ 
found as 



the nonlinear transmission spectrum in the vicinity of the resonance can be 



Tioo) 



Pr 



"W ^ 1^0 _|_ -Pout 
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-1 



(71) 



For a fixed value of the frequency detuning, the formula ( 71 ) can be used to compute the 



hysteresis curve. It can be checked that this curve is equivalent to the polynomial of the third 
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degree with real coefficients, and thus it can describe only bistable resonances. More complex 
structures which demonstrate multistable behavior can be constructed by combining several 



strictly bistable microcavities 11 . To treat the nonlinear properties of such structures, it is 



useful to define the nonlinear transfer matrix T,, 





CO — CUq \v+ ± v\ 

^ + ^^ — I 

7 ^0 



1 



±1 
-1 



(72) 



(73) 



where I is the identity matrix. The transfer matrices of single microcavities can be multiplied 
producing the total transfer matrix of the structure. This leads to a hysteresis curve of more 
complicated shape which can be different for the opposite directions of incidence because 
the transfer matrices do not commute and the order of multiplication plays an important 
role IgI. 



5. Numerical examples 

5. A. Bragg gratings with symmetrically placed defect 

Multilayered structures can be considered as a one-dimensional (ID) realization of on-channel 
microcavities and are particularly suitable to check the accuracy of the CMT equations. As a 
test case, we use a Bragg structure with a symmetrically placed defect which can be described 
by the symbolic formula (HL)p(LH)^. In what follows, the letters 'L' and 'H' correspond 
to quarter-wave layers of polydiacetylene 9-BCMU with linear (nonlinear) refractive index 
nL = 1.55 (nsL = 2.5 x 10"^ cmVMW) and rutile with rin = 2.3 (nan = 10"^ cmVMW), 
respectively [23, 



24 



. The quarter wave condition is set to Aq = 0.7 /im 
nhdh = niidn = Aq/4, 



(74) 



which gives the thicknesses of the layers di, = 112 nm and du = 76 nm. The main advantage 
of this structure is that it has a well-defined resonance at Uq = 27rc/Aq which is surrounded 
by band gap regions. Due to the mirror symmetry of the structure, this resonance always 
shows perfect transmission, and its half-width can be adjusted by the number of periods p 
in the Bragg mirrors. 



The formulas for the decay rate (31 ) and the characteristic intensity (66 ) of the microcavity 
at the resonance Aq = 2ttc/uq can be simplified in the ID case to 

1 -1 



7 



'0 
-L 



2vr/Ao) en2EQdx 



1 -1 



(75) 
(76) 
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Instead of the decay rate 7, it is often convenient to use a dimensionless quality factor defined 
as Q = uo/ (27) 

Q = (tt/Ao) {"^eE^dx. (77) 



In some simple cases, the integration can be performed analytically, and an explicit formula 
for the quality factor can be obtained [see Appendix ]. 



For the structure (IIL)^(LH)^, the quality factor computed by Eq. (77) is Q = 2060 and the 



characteristic intensity according to Eq. (76) is Jq = 0.05695 MW/cm . The two parameters 
together with Eq. ( 71 ) fully determine the hysteresis and the nonlinear transmission spectrum 
of the structure in the vicinity of the resonant frequency uq = Uq. The quality factor can be 
also computed with the linear transfer matrix by finding the full-width at half-maximum of 
the resonance cjq/ Acjfwhm = 2058. Therefore, the accuracy of the CMT equations in this 
case can be estimated as 0.1%. 

It is worth noting that the CMT parameters of a similar structure (LH)^(HL)^ are different. 
It has a smaller quality factor Q = 581.3, and a significantly larger characteristic intensity 
Jo = 1.593 MW/cm^. This can be explained by a different localization of the electric field 
in the structure [Fig. [3]^a,b)]. Nevertheless, the usage of normalized units ensures that both 
structures have the same shape of the hysteresis and the nonlinear transmission spectrum 
[Fig-ic-f)]. 

The Maxwell equations for nonlinear multilayered structures 

d^Ey = i{u/c)H,, (78) 
d^H, = i{u/c)e{l + 2n2\Ey\^)Ey (79) 



can be also solved by using strictly numerical methods 25 , 26 or semi-analytical techniques 



27,28. A detailed comparison with the results obtained by the Runge-Kutta method is 



presented in Fig. |3[c-f). 

5.B. Thue-Morse multilayered structures 

As a more complex example, we consider a Thue-Morse quasicrystal. It has a nonperiodic 
arrangement of layers which is governed by a deterministic set of inflation rules and features a 
number of pseudo band gap regions with resonances of complete transmission flT] . We choose 
one of such resonances which is located at uq = 0.705465 Uq. It has two localization centers in 
the field profile so that the full structure can be divided into two parts which can be treated 
as coupled microcavities [Fig. |4]^a)]. These microcavities, which will be denoted as a and (3, 
have the same resonant frequencies Ua = oJ/s = ujq, but their parity and CMT parameters are 
different. For the left part, the quality factor is Qa = 318.2 and the characteristic intensity 
is la = 5.147 MW/cm^, while for the right part Q/^ = 1110 and J/? = 0.4078 MW/cm^. 
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The nonlinear response of these microcavities is quahtatively similar and can be described 
by the same analytic formula (71) [Fig. |4]^b,c)]. The discrepancy with the numerical results 
is noticeable only for the microcavity a which has a relatively small quality factor. Since 
the accuracy is worse on the higher frequency side of the resonance where the band gap 
is less pronounced, this suggests that the influence of other resonances causes additional 
perturbations. 

The coupling between microcavities leads to more complex hysteresis curves and nonlinear 
transmission spectra [Fig.|4|^d,e)]. They can be obtained by multiplying the nonlinear transfer 
matrices of the microcavities (73) for a fixed output intensity and then restoring the input 
intensity. The order of multiplication plays an important role in the nonlinear case because 
the transfer matrices do not commute and the result strongly depends on the direction of 
incidence. Apart from the nonreciprocal behavior, there is also the possibility of multistable 
behavior since the hysteresis curve in the case of coupled microcavities is described by a 
polynomial of a higher degree. It is very important that a simple model based on the CMT 
equations is able not only to explain the nonlinear properties of complex resonances like this 
one, but also shows a good quantitative agreement with computationally intensive numerical 
methods. 



6. Conclusions 

The phenomenological CMT is a very efficient tool for studying the nonlinear behavior of 
microcavities both in the frequency and time domain. It considers the interaction between 
microcavity and waveguide modes in a way that is similar to the scattering formalism. 
Therefore, the complex wave dynamics can be separated from a relatively simple picture of 
coupling, and this gives a significant advantage in comparison to strictly numerical methods. 
The dynamical properties of the microcavities can be fully determined by a small set of 
parameters which includes the decay rate, coupling coefficients and characteristic intensities. 

By using on-channel microcavities with two coupling ports as an example, we provided 
for the first time a systematic derivation of the CMT equations starting directly from the 
Maxwell equations and obtained the explicit formulas for all phenomenological parameters. 
Our derivation is particularly suitable for microcavities embedded in photonic crystal waveg- 
uides of various dimensionality and multilayered structures. The accuracy of the results de- 
pends on the quality factor of a specific resonance and is mostly limited by the infiuence of 
other resonances. 
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Appendix A: Quality factors of Bragg gratings with symmetrically placed defects 

By using the fact that the energy density is a constant in each layer of the structure, the 



formula for the quality factor (77) can be rewritten as 

/An — 



^0) 



where dk is the thickness of the layer k with the refractive index ra^, and the sum is taken 
over all layers in the structure. Fields on opposite sides of the layer k can be related by the 
characteristic matrix 





Mi 



irij^ ^ sin 



(82) 



COS^fc 

^ mfcsin^fc cos^fc 

UkdkUj/c, or ^k = (7i"w)/(2u;q) if all layers satisfy the quarter wave condition 
(74) at the frequency Uq. The M- matrix for the single period of the Bragg gratings in the 
structure (HL)p(LH)p can be obtained as a multiplication of M-matrices corresponding to 
layers 'L' and 'H'. It takes a particularly simple form at the resonance uq 



where 



MlMh 



3) 



Wo 

nn/riL 
ni,/nu 

which shows that the fields are exponentially growing or decaying towards the center of 
the structure as (wh/'t-l)^, where p is the number of periods in the Bragg mirrors. The 
contribution of each period to the quality factor is 



TT 



\Ek\' + 



which makes in total 



P-I r / N 2k+l , 

Q=>H + nOE {-) +i 



2k 



^4) 



^5) 



The sum of the geometric progressions can be found as ^ 



k=0 ' 



(1 - rP)/{l - r) and 



keeping only the largest term leads to the following formula for the quality factor of the 
structure (HL)p(LH)p 



Q 



4(nH - nL) 

It is worth noting that the quality factor of a similar structure (LH)p(HL)p 



Q 



TT 



4(nH - rii,) V^L/ 



(86) 



^7) 



is smaller in tih^^l times. 
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Fig. 1. (Color online), (a) Photonic crystal waveguide with an embedded 
defect which behaves as a microcavity. The region where the electromagnetic 
filed is negligible has a dark gray background. The integration volume for 
the eigenvalue problem is shown by the dashed line. The dominant flow of 
energy goes through the reference planes at x = and x = L. (b) A sketch 
of the model used in the scattering formalism. The amplitudes of the ingoing 
waves («+ and v^) are related to the outgoing waves («_ and f_) through the 
amplitude of the microcavity mode {A). 
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Fig. 2. (Color online). An example of two fundamental solutions that exist at 
the same resonance frequency and satisfy the boundary conditions of standing 
waves with zero energy flow. The electric and magnetic fields are shown by 
solid and dashed hues, respectively. For one of the solutions they are exponen- 
tially growing towards the defect region in the middle (a) while for the other 
independent solution they are exponentially decaying (b). Background shows 
alternation of layers with a higher ('H') and lower ('L') index of refraction 
inside the structure. 
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Input intensity, 1 1 Iq Frequency detuning, Ao) / y 




Input intensity, 1 1 Iq Frequency detuning, Ao) / y 



Fig. 3. (Color online). The structures (HL)^(LH)8 (a) and (LH)S(HL)^ (b) 
show a strong resonance at the quarter wave frequency Uq = Uq, but have a 
different distribution of the electric field. The hysteresis of transmission for 
a fixed value of the frequency detuning (c, e) can be computed by Eq. (71), 
which follows from the CMT (solid lines), and by solving the Maxwell equa- 
tions (78, 79) directly with the Runge-Kutta method (circles). Not all parts 
of the hysteresis curves are stable, and arrows show how the switching be- 
tween different branches occurs. The usage of normalized units ensures that 
both structures have the same shape of the hysteresis curve. The results of 
the two methods are in good agreement (the comparison for (HL)'^(LH)^ is 
shown, (LH)'^(HL)^ is similar). The same methods can be applied to compute 
the nonlinear transmission spectrum for a fixed value of the input intensity 
(d, f). 
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Fig. 4. (Color online). An example of the resonance (the Thue-Morse struc- 
ture of the 7th generation number at the frequency ojq = 0.705465 Wq) which 
has two localization centers in the electric field profile (a) denoted as a and /3. 
The hysteresis (b) and nonlinear transmission spectrum (c) of these parts can 
be computed by Eq. (73), which follows from the CMT (solid lines), and by 
solving the Maxwell equations ( 78 , 79 ) directly with the Runge-Kutta method 
(circles). The discrepancy between the two methods is noticeable only for the 
part a, because it has a relatively small localization strength. The Thue-Morse 
structure, which combines the parts a and /3 together, shows the nonrecip- 
rocal behavior. For the same frequency u = 0.7046 Uq and input intensity 
J = 1.5 MW/cm^ as used in (b) and (c), respectively, the hysteresis (c) and 
nonlinear transmission spectrum (e) depend on the propagation direction (a/? 
or f3a). 
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